## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Introduction

As the climate changes, spring phenological events like budburst and flowering will advance, especially for plants active in rapid seasonal transitions and short growing seasons [@Pau2011], like many high elevation and latitude conifers. This effect is already obvious in many species [@Parmesan2006; @Franks2014]. Changes in pollination phenology can affect fecundity, gene flow, and even range size in a species and have effects on dependent species [@Inouye2008; @Chuine2001a].

Conifers are a big part of the enormous northern northern hemisphere forests and they have wide ranges with lots of local adaptation. Common garden experiments and genetic work reveal extensive local adaptation in many forest tree species, especially boreal and temperate conifers (reviewed in @Alberto2013a). A locally adapted population only grows optimally in a subset of the range and may tolerate a more limited climatic range than the species as a whole. In northern hemisphere conifers, local adaptation often reflects strong trade-offs between avoidance of cold damage and competitive height growth (summarized in @Aitken2008a).

Coniferous forest trees are wind pollinatated wtih pollination possible over large distances Pollen is shed from male strobili and must arrive at receptive female strobili for successful pollination. Shifts in the timing of pollen shed and cone receptivity (pollination phenology) in conifers could lead to gene flow changes that hinder or promote adaptation under climate change, decrease fitness, and even affect reforestation via seed production declines. [Also affects gene flow now and is important for understanding current spatial genetic structure and local adaptation]

Lodgepole pine is a good representative of these issues - an economically and ecologically important tree species facing multiple threats from climate change [@Schneider2010; @Sambaraju2012; @Hamann2006]. Lodgepole pine has a very large geographical distribution (across 33\(^\circ\) latitude and 31\(^\circ\) longitude) encompassing a wide range of climates and soils (Figure ) with widespread and significant local adaptation in many traits. For example, populations from both northern interior British Columbia and northern Idaho can survive in areas with mean annual temperatures between -4 and 6 \(^\circ\mathrm{C}\), but the northern British Columbia population survives best where mean annual temperatures are ~ 1 \(^\circ\mathrm{C}\) and the Idaho population best at ~ 4 \(^\circ\mathrm{C}\) [@Rehfeldt1999a]. Local adaptation in lodgepole pine can be observed even at relatively small spatial scales when topographic variability is high: in a reciprocal transplant experiment, growth declines were observed when moving high elevation populations just 100m in elevation [@Rehfeldt1983]. [Summarize/cite adaptree work]

Rangemap of lodgepole pine. After @Little1971 .

Rangemap of lodgepole pine. After @Little1971 .

Pollen is an important vector for identifying gene flow in lodgepole pine because outcrossing is common, pollen dispersal is extensive, and seed dispersal is relatively limited [@Ennos1994]. There is evidence for spatially varying levels of gene flow in the species as populations from areas with higher regional climate heterogeneity have higher genetic variance [@Yeaman2006]. Pollination phenology could control this.

Pollination phenology determines which populations can exchange genes, but predicting the timing of pollen shed and cone receptivity has not been done in lodgepole pine. Pollination phenology examples are uncommon. [Talk about Sarvas investigations for mechanistic background, why using temperature, etc.] Simple heat accumulation thresholds (Pinus taeda) [@Boyer1978] or elevation (Pinus flexilis) [@Schuster1989] were used previously to explain or predict pollen shed in limited spacial and temporal contexts. @Owens2006 reports that lodgepole pine pollen shed and cone receptivity occur when degree days reach about 500 at a threshold of 5 \(^\circ\mathrm{C}\), but this is the only report of pollination phenology modeling for lodgepole pine and limited details are provided. Models of lodgepole pine vegetative phenology, on the other hand, are better represented in the literature (e.g. @Chuine2001), and pollen shed and cone receptivity are not expected to have additional or more complex triggers or model forms than budburst [@Chuine2013a].

Predicting pollination phenology will also have practical benefits. Seed orchard managers in British Columbia are particularly concerned about protandry, when all pollen in an area is shed before cones become receptive (personal communication, Chris Walsh, Former Seed Orchard Manager, Kalamalka Seed Orchards, February 13, 2013). Protandry occurs in particularly hot and dry years [@Owens2005]. If this pattern holds, protandry could become more common in natural populations, leaving some populations pollen limited and likely hampering local adaption. [Outcrossing! Inbreeding!]

Crucial to the scale of this project, ~15 years of pollination phenology data are available from seed orchard sites producing seed for reforestation for lodgepole pine. Seed orchards in British Columbia produce large amounts of tree seed for reforestation from parent trees sourced from provenances around the province. Interior lodgepole pine (Pinus contorta var latifolia) is planted more than any other species at rates of [XXX] million/year as of [201X]. In the late 90s, it became apparent that seed yields (filled seeds per cone) at seed orchards in the northern Okanagan were much lower than in orchards near Prince George. To plan for future seed production and orchard establishment and management, studies were undertaken to understand why seed yields at north Okanagan orchards were lower than at the Prince George Tree Improvement Station orchards.Pollination phenology data was collected at the Prince George Tree Improvement Station in British Columbia beginning in 1997 and collection at many other BC tree orchard sites began in 2006 under the Forest Genetics Council’s Operational Tree Improvement Program 0722 [@Webber2007] . Pollination phenology data recorded as part of this program included pollen concentrations and dates of pollen shed and cone receptivity.

I will use pollination phenology and temperature data to fit a mechanistic model predicting pollen shed and cone receptivity across the entire range of Pinus contorta var. latifolia and consider how the date and length of pollination phenology periods vary normally and under several climate change scenarios. Specifically, I will answer 1) What is the relationship between temperature and pollen shed and temperature and cone receptivity timing and length? 2) Does provenance affect that relationship? 3) How many days will pollination phenology shift under the next 30, 60, and 90 years of climate change 4) Will protandry become more common?

Methods

Study species

While lodgepole pine has four subspecies, this work concerns only Pinus contorta subsp. latifolia.

Data

Phenology Data

My aim was to model pollination phenology in lodgepole pine and calculate synchrony across the distribution. To determine the timing of pollen shed and cone receptivity, I modeled the forcing requirements for pollination phenology in lodgepole pine and investigated differences in forcing requirements between males and females and among provenances. I then applied the model to weather data across the lodgepole pine distribution to determine the distance between populations with synchronous pollen shed and cone receptivity.

Sites

Trees selected from across the British Columbia portion of the lodgepole pine range are grown in both large scale provenance trials and seed orchards, which are similar to common gardens, as part of tree breeding and seed production programs. Between 1997 and 2011, flowering phenology of lodgepole pine was recorded at 8 seed orchard sites in the interior of British Columbia. I contacted Seed Orchard Managers and other forestry professionals across British Columbia in 2012 and received pollination phenology data from C. Walsh, previously at Kalamalka Seed Orchards (now retired), R. Wagner at the Prince George Tree Improvement Station, and J.E. Webber previously at the Glyn Road Research Station (now retired). 4 of the sites are clustered near to one another, but sites span about 5 degrees of latitude and are at elevations from [X] to [X].

Map of seed orchard locations. Boxed area in map on left is shown in greater detail on the right.

Map of seed orchard locations. Boxed area in map on left is shown in greater detail on the right.

Provenances

Trees grown at the seed orchard sites were sourced from 6 biogeoclimate regions known as Seed Planning Units. Trees with the same SPU provenance are grown together in an orchard at a given site. Genotypes (labelled with a Clone number in the data) in the orchards are represented by multiple ramets.

I obtained records for 17 of the 26 lodgepole pine seed orchards in British Columbia. Orchards include the offspring of trees from 6 of the 7 BC seed planning zones (SPZs) (Figure ) grown at at 7 sites across BC. SPZs are biogeoclimatic and political units used for seed planting purposes by British Columbia. SPZs are divided into elevation bands called Seed Planning Units (SPUs), which form this project’s provenances. Thirteen out of 17 orchards in my data set are first generation orchards and should faithfully represent their provenances. These first generation orchards represent 6 provenances at 5 sites. Second generation orchards have been selectively crossed and this may skew the mean or variance of phenology for a provenance if pollination phenology varies by provenance.

Most provenances are represented at 2 to 3 sites and have at least three years of data at a given site spanning 1997-2012 (Figure ). The Prince George Tree Improvement Station (PGTIS) provides a continuous 15-year record of its three orchards’ phenology.

Contingency table of years of data for Seed Planning Units (rows) and Seed Orchard Sites (columns). Seed Planning Zones, used as provenances in this project, are usually represented at multiple years and multiple sites. There is particularly good representation at PGTIS.

Contingency table of years of data for Seed Planning Units (rows) and Seed Orchard Sites (columns). Seed Planning Zones, used as provenances in this project, are usually represented at multiple years and multiple sites. There is particularly good representation at PGTIS.

Map of Seed Planning Units (SPUs). Seed planning units are biogeoclimatic and political units used for seed planting purposes by British Columbia. Seed planning units form this project’s provenances. High, Low, and Mid refer to elevational bands. Data is also available for East Kootenay Low, but will likely not be included in any analysis as it includes only one year at one site.

Map of Seed Planning Units (SPUs). Seed planning units are biogeoclimatic and political units used for seed planting purposes by British Columbia. Seed planning units form this project’s provenances. High, Low, and Mid refer to elevational bands. Data is also available for East Kootenay Low, but will likely not be included in any analysis as it includes only one year at one site.

[table with Site Columns and SPU rows with years of data as values][table of number of trees/clones in a given year for a given site/spu combo]

Phenology scoring protocol

Protocol C in [@Woods1996] was used as the basis for collecting pollen shed and cone receptivity data, though operational constraints led to some modifications. Workers monitored seed orchards for the beginning of pollen shed and cone receptivity. ~15 clones, usually represented by 2 trees each, were selected for specific observations. When the active period seemed to be starting, workers went into the orchard every few days to make observations on the selected trees.

Stages of pollen and seed cone development are described by [Owens & Molder 1984 and updated in Owens2006] and were used as a general guide for determining the phenological state of pollen and seed cones. Pollen cones are flowering when tapping causes pollen to be released and seed cones are flowering when there are gaps between the bract-scale complexes. Pollination drops are also produced during flowering, though they recede midday if pollen is not present [@Owens2006].

“Flowering” states were recorded for both pollen shed and cone receptivity at the level of each tree. Protocol C recommends marking the dates when 20% of the cones on a tree have begun flowering and when 80% of the cones on a tree have finished flowering. Operationally, there was some subjectivity and tree-level states for each cone type should be interpreted as “starting flowering” and “finished flowering.” There is some subjectivity here.

[DESCRIBE SURVEY PERIODS, MAYBE MAKE A TABLE OR GRAPH OF OBSERVATION DATES?]

Observations were not made every day and survey periods varied in length. At Prince George, not all trees have complete phenological records, e.g if a tree is not flowering on the first day of observation, the start date is unknown.

Harmonization and cleaning

Describe data transcription and cleaning.

There were some differences in how data was recorded at the Prince George Tree Improvement Station versus the other sites. At Prince George, trees were marked as flowering or not flowering on each day of observation. At other sites, only the first day observed flowering and the first day observed finished flowering were recorded. I inferred that trees were not yet flowering on observation days prior to their first recorded flowering date. I cleaned and harmonized the data for analysis in a single model using R scripts provided so that, for pollen cones and seed cones, stage 1 = not yet flowering, stage 2 = flowering, and stage 3 finished flowering.

There are three phenological stages of interest each for a tree’s male and female cones

    1. The cones have not yet flowered
    1. The cones are flowering
    1. The cones have finished flowering

Phenophases in the field were recorded using different symbol sets and resolutions. I assigned each symbol to one of the phenophases above. Trees that did not produce cones are assigned phenological stage 0.

Phenophase Symbols Male.Cones Female.Cones
0 0 none produced none produced
1 1, 2.5, - not yet shedding not yet receptive
2 3, 3.5, 4, 4.5, 5, pollenshed20, receptive20 shedding receptive
3 -, receptive80, pollenshed80 finished shedding no longer receptive

Weather data

Daily weather data at seed orchard sites was extracted from PNWNAmet, a daily gridded meteorological dataset at 1/16 \(\circ\) over northwest North America [@pacificclimateimpactsconsortiumPNWNAmet194520122014]. The closest point in the PNWNAmet was used for each station. Mean daily temperature was calculated from the minimum and maximum daily temperatures.

Map of Sites with nearest climate gridpoints

Map of Sites with nearest climate gridpoints

Forcing units

The relative effect of mean daily temperature \(t\) on the rate of forcing was calculated with

\[R(t) = \frac{28.4}{1+ e^{-0.185*t-18.4}}\] which describes the relationship between temperature and development rate. This equation for forcing units is based on experimental work from Risto Sarvas with in temperate forest tree species [@Sarvas1972, @hanninenModellingBudDormancy1990]. @Chuine1999 found this calculation of forcing units superior to degree day calculations.
Forcing Units on day \(d\) (\(f_d\)) are the sum of the relative temperature effect from January 1 (ordinal date \(1\)) to day \(d\).

\[f(d) = \sum_{i=1}^d R_d(x)\]

The rate of forcing was experimentally determined by Sarvas [@Sarvas1972] and verified to perform better than growing degree day models by Chuine [@Chuine1999]

Phenology Modeling

Actual methods

I modeled the derived phenological states with a Bayesian multilevel ordered logistic model fit in Stan with the no-u-turn sampler [@standevelopmentteamRStanInterfaceStan2018]. An ordinal logistic model

For phenological states \(k \in \{1,2,3\}\), the probability of

The ordered logistic probability mass function is \[\text{OrderedLogistic}(k~|~\eta,c) = \left\{ \begin{array}{ll} 1 - \text{logistic}(\eta - c_1) & \text{if } k = 1, \\[4pt] \text{logistic}(\eta - c_1) - \text{logistic}(\eta - c_2) & \text{if } k=2, \text{and} \\[4pt] \text{logistic}(\eta - c_3) - 0 & \text{if } k = 3. \end{array} \right.\]

\[\begin{array}{rlr} S_i & \sim \mathrm{OrderedLogistic}(\eta,\kappa) & \text{probability of data}\\ \eta_i & = (\beta_{site} +\beta_{prov} + \beta_{clone} + \beta_{year}) f_i & \text{linear model}\\ \kappa_k & \sim \mathrm{gamma}(7.5,1) & \text{fixed priors}\\ \beta & \sim \mathrm{exponential}(2) & \text{}\\ \beta_{site} & \sim \mathrm{normal}(0,\sigma_{site}) & \text{adaptive priors}\\ \beta_{prov} & \sim \mathrm{normal}(0,\sigma_{prov}) & \text{}\\ \beta_{year} & \sim \mathrm{normal}(0, \sigma_{year}) \\ \beta_{clone} & \sim \mathrm{normal}(0, \sigma_{clone}) \\ \sigma_{site} & \sim \mathrm{exponential}(3) & \text{hyperpriors}\\ \sigma_{prov} & \sim \mathrm{exponential}(3) \\ \sigma_{year} & \sim \mathrm{exponential}(3) \\ \sigma_{clone} & \sim \mathrm{exponential}(3) \\ \end{array}\]

The point at which a tree is 50% likely to have transitioned from one phenophase to another can be calculated by dividing the cutpoint for the transition by the slope.

\[\begin{array}{ll} \bf{H}_1 &= \kappa_1/(\beta_{site} +\beta_{prov} + \beta_{clone} + \beta_{year})\\ \bf{H}_2 &= \kappa_2/(\beta_{site} +\beta_{prov} + \beta_{clone} + \beta_{year})\\ \end{array}\]

Another way to think about the probabilities in this model:

\[\begin{array}{rll} S_i&\sim \mathrm{Categorical}(\bf{p}) &\text{probability of data} \\ p_1 & = q_1 &\text{probabilities of each value of } k\\ p_2 & = q_2 - q_1 \\ p_3 & = 1-q_3 \\ logit(q_k) &= \kappa_k - \eta_i &\text{cumulative logit link} \\ \eta_i&= (\beta_{site} +\beta_{prov} + \beta_{clone} + \beta_{year}) f_i & \text{linear model} &\text{linear model}\\ \\ \end{array}\]

To quantify accuracy of the model, I will calculate the difference in variance within and between populations. Only phenological differences larger than the within population variance will be considered. The maps and parameters of these models can be compared to answer the following questions

Differences between provenances

Were calculated according to [@vandermijnsbruggeVariationBudBurst2015] - “difference in days in which half the plants has reached the same score level”

\[\mathrm{H50}_{prov_i} - \mathrm{H50}_{prov_j} = \beta_{prov_j} - \beta_{prov_2}/\beta\]

Assumptions

Chilling requirements are met and chilling and forcing periods do not overlap. Forcing is a function of temperature only, not photoperiod or chilling. Both phases occur at the same rate.

Results

  1. How strong are the genetic vs. environmental effects on pollination phenology?

To answer this question I will compare the parameters of models fit by provenance. If different provenances have significantly different model parameters - that is, have different heat sum requirements for pollen shed and cone receptivity by provenance, then the rest of my project will be limited to the provenances I have data from. If not, I can use the same models across the entire range when predicting phenological events.

The budburst and spring pollination phenology of temperate tree species is determined by temperature, and there may be some regional variation in the temperatures required for pollen shed and cone receptivity; in lodgepole pine, the threshold for shoot elongation was 5.1 \(^\circ\mathrm{C}\) for coastal and interior provenances, but 6.5 \(^\circ\mathrm{C}\) for northern provenances [@Chuine2001].

However, seedling spring vegetative phenology data from four growth chamber temperature and moisture treatments in the AdapTree project suggests that budburst shows substantial clinal variation (> 10% of phenotypic variance explained) only in the coldest treatment tested (MAT of 1 \(^\circ\mathrm{C}\)), and much colder than the seed orchard locations (Figure , Figure , @Liepe2014 (submitted)). Analysis of phenotypic data from an outdoor common garden in Vancouver (personal commmunication, Ian Maclachlan, AdapTree Graduate Student, January, 1 2015) also shows that very little variation exists among provenances for spring phenology, regardless of whether seedlings originate from wild-stand or selectively bred seed orchard seedlots.

Variance in budbreak and budset for lodgepole pine explained by provenance climate. Seedlings from 281 locations across lodgepole’s BC range were grown in growth chambers under four treatments: approximated seasonal and daily variation in temperature for geographic locations with a MAT of 1, 6, and 11 ^\circ\mathrm{C}, with the warmest treatment having a dry and wet version. Pearson’s correlation coefficient between phenotypic traits in lodgepole pine and provenance climates are presented. The vertical line represents the critical r^2 value after an adjustment for multiple inferences. From @Liepe2014 (submitted).

Variance in budbreak and budset for lodgepole pine explained by provenance climate. Seedlings from 281 locations across lodgepole’s BC range were grown in growth chambers under four treatments: approximated seasonal and daily variation in temperature for geographic locations with a MAT of 1, 6, and 11 \(^\circ\mathrm{C}\), with the warmest treatment having a dry and wet version. Pearson’s correlation coefficient between phenotypic traits in lodgepole pine and provenance climates are presented. The vertical line represents the critical \(r^2\) value after an adjustment for multiple inferences. From @Liepe2014 (submitted).

  1. How does faster spring warm-up like we expect under climate change affect within population mating success and frequency of outcrossing?

I expect that faster spring warm-up will condense as well as advance growth initiation dates and also shorten pollination phenological periods by speeding development prior to and during the phenological period. This should further synchronize populations in similar climates and act as a barrier to breeding with populations in more different climates. On the other hand, it could favor outcrossing over selfing by causing protandry - pollen shed in advance of cone receptivity. Based on anecdotal reports from seed orchards managers, I expect protandry (and thus outcrossing) to increase under climate change.

I will examine the variability of the start date across the range (\(t_1\)) to see if it shrinks in warmer years and as the climate changes. I will also calculate overlap between within population cone receptivity and pollen shed under current conditions and with climate warming and see if the heat sum requirements for cone receptivity and pollen shed are different.

Discussion

Grafting bias?

Problems from breeding

There are two types of orchards: first generation and advanced. Each orchard contains trees that are the descendants of seeds and scions collected from mother trees in one particular provenance. First generation orchards have not undergone selective breeding but may be subject to selection through the process of choosing mother trees (called “plus” trees by tree breeders), testing and selection of mother tree offspring, or culling of unfavorable trees in the orchard. First generation orchard trees are clones of mother trees and their offspring grafted onto rootstock. There are three types of first generation orchards subject to different levels of selection: 1.0, 1.5, and 1.75 generation orchards. All orchards are subject to selection based on mother tree selection. In 1.5 generation orchards, superior offspring of mother trees are selected for the orchard or poor trees are culled from the orchard. This selection may be strong; more than 50% of offspring can be rejected in this process due to poor growth form [@Ukrainetz2011]. In 1.75 generation orchards, data from 10 year tests of mother tree offspring are used to select 40 of the best genotypes for a given provenance. Offspring from controlled crosses between trees in first generation orchards are tested in field trials and the best trees are then cloned and grafted onto root stock in advanced generation orchards.

orchards_in_data <- read.csv("~/Documents/research_phd/data/PhenologyAndPollenCounts/data_formatted_and_derived/derived_phenophase.csv", header=TRUE, stringsAsFactors = FALSE) %>% 
    select(Orchard, Year, ) %>%
    group_by(Orchard) %>%
    dplyr::summarise(Years = n_distinct(Year))
orchard_gen <- read.csv("~/Documents/research_phd/data/OrchardInfo/OrchardGen.csv", header=TRUE, stringsAsFactors = FALSE)

orchard_meta <- dplyr::left_join(orchards_in_data, orchard_gen) %>%
    unique()
## Joining, by = "Orchard"
knitr::kable(orchard_meta)
Orchard Years Generation
218 2 1.5
219 1 1.5
220 15 1.75
222 2 1.5
223 14 1.75
228 15 1.75
230 4 1.5
234 2 Advanced
237 3 Advanced
240 2 Unknown
307 4 1.75
308 1 1.75
310 1 1.5
311 1 1.5
313 1 1.5
338 2 Advanced
339 2 1

selection only from a portion of the range

Provenances and sites included in my data exclude the coldest and wettest parts of the range. Figure shows the 1961-1990 climate normal mean annual temperature (MAT) and mean annual precipitation (MAP) for data from all lodgepole pine provenances in BC, provenances included in my data set, and sites where provenances were grown. Orchard sites are generally much warmer and drier than many locations where lodgepole pine grow. I may have difficulty extrapolating into the coolest and wettest parts of the lodgepole pine range, depending on actual yearly conditions at sites for which I have data. However, this may be advantageous when it comes to predicting phenology under climate change.

Provenance (SPU) climates of data included in this proposal, provenances not included in this proposal, and sites where provenances were grown. MAT is mean annual temperature and MAP is mean annual precipitation.

Provenance (SPU) climates of data included in this proposal, provenances not included in this proposal, and sites where provenances were grown. MAT is mean annual temperature and MAP is mean annual precipitation.

no chilling

Previous models for lodgepole pine have fitted spring growth phenology models without considering chilling [@Chuine2001] and, in general, budburst models for boreal tree species have not required chilling to give accurate predictions [@Linkosalo2000]. For lodgepole pine growing in natural environments, chilling is always met and should continue to be met under the next century of climate change; in the AdapTree project, lodgepole pine seedlings grown for three to four weeks at 4 \(^\circ\mathrm{C}\) met their chilling requirements. Thus, I expect to be able to disregard the state of chilling and fit only Equations and ). Once the model is parameterized I will be able to calculate the timing of pollen shed and cone receptivity across the lodgepole pine species range using location specific climate data from seed orchards and ClimateWNA [@Wang2012]. I will validate the model by comparing predictions to phenological measurements from the National Phenology Network.